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Abstract: For the most popular sequential change detection rules such as CUSUM, 
EWMA, and the Shiryaev-Roberts test, we develop integral equations and a concise 
numerical method to compute a number of performance metrics, including average 
detection delay and average time to false alarm. We pay special attention to the 
Shiryaev-Roberts procedure and evaluate its performance for various initialization 
strategies. Regarding the randomized initialization variant proposed by PoUak, 
known to be asymptotically optimal of order-3, we offer a means for numerically 
computing the quasi-stationary distribution of the Shiryaev-Roberts statistic that 
is the distribution of the initializing random variable, thus making this test appli- 
cable in practice. A significant side-product of our computational technique is the 
observation that deterministic initializations of the Shiryaev-Roberts procedure can 
also enjoy the same order-3 optimality property as PoUak's randomized test and, 
after careful selection, even uniformly outperform it. 
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1. Introduction 

Change-point problems deal with anomaly detection or more generally detec- 
tion of changes in the statistical behavior of processes. This problem has an 
enormous spectrum of important applications, including biomedical signal and 
image processing, quality control engineering, financial markets, link failure de- 
tection in communication networks, intrusion detection in computer networks 
and security systems, detection and tracking of covert hostile activities, chem- 
ical or biological warfare agent detection systems (as a protection tool against 
terrorist attacks), detection of the onset of an epidemic, failure detection in man- 
ufacturing systems and large machines, target detection in surveillance systems, 
econometrics, seismology, navigation, speech segmentation, and the analysis of 
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historical texts. See, e.g., Willsky (1976), Basseville and Nikiforov (1993), Mac- 
Neill (1993), Baron (2002), Galstyan, Mitra, and Cohen (2007), Kent (2000), 
Tartakovsky (1991), Tartakovsky and Ivanova (1992), Tartakovsky and VeeravaUi 
(2004), Tartakovsky, Li, and Yaralov (2003), Tartakovsky, Rozovskii, Blazek, and 
Kim (2006), Wang, Zhang, and Shin (2002). In all of these applications, sensors 
monitoring the environment take observations that undergo a change in distri- 
bution in response to changes and anomalies in the environment or changes in 
patterns of a certain behavior. The observations are obtained sequentially and, 
as long as their behavior is consistent with the normal state, one is content to 
let the process continue. If the state changes, then one is interested in detecting 
the change as soon as possible while minimizing false detections. 

Let {Xn}n>i denote observations that arc obtained sequentially, and let F^o 
and Po be the probability measures before and after the change. Let Fj- and E,- 
denote the probability measure and the expectation induced when the time of 
change is r > 0. We use the convention that r is the last time instant where 
the observations follow the nominal (pre-change) regime, so the first observation 
under the alternative measure is at time r + 1. Thus Pq means that the change 
took place before any observations were taken and all observations are under the 
alternative regime, whereas Pqo stands for the scenario in which the change is at 
infinity (i.e., does not occur) and all observations are under the nominal regime. 

A sequential change detection procedure is identified with a stopping time 
T that is adapted to the filtration {J-'n}n>o, where J^o is the trivial cr-algebra 
and, foi n > 1, = cr{Xi, . . . ,Xn} is the cr-algebra generated by the first n 
observations. Thus the event {T < n} belongs to J^n- Since the change occurs 
at an unknown instant, the objective is to detect it as quickly as possible while 
avoiding frequent false alarms. Therefore, the design of the quickest change-point 
detection procedure involves optimizing the trade-off between two kinds of per- 
formance measures, one being a measure of detection delay and the other being 
a measure of the frequency of false alarms. Regarding the latter, the false alarm 
rate is usually measured by the average time to false alarm EoofT"], commonly 
referred to as the average run length (ARL) to false alarm. For detection delay, 
two major non-Bayesian criteria have been proposed in the literature. The first, 
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due to Lorden (1971), is 

Ji^{T) =supesssupE^[(T-r)+|J^^], (1.1) 

r>0 

where = max{3;,0}; and the second, due to Pollak (1985), is 

Jp(r) =supE,[r-r|r > r], (1.2) 

T>0 

where Mt-[T — t\T > r] is the (conditional) average delay to detection for the fixed 
point of change < r < oo. As discussed in Moustakides (2008), Lorden's perfor- 
mance measure is appropriate for problems in which the change-point mechanism 
takes into account the observations for deciding about imposing the change. Pol- 
lak's measure, on the other hand, assumes that the change is imposed by a source 
independent from the observations. Both cases are equally important, referring 
to completely different classes of change-point applications. 

Whether we use Lorden's or Pollak's measure, in finding an optimal stopping 
time one would normally be interested in minimizing Jh{T) or J-p{T) and maxi- 
mizing, at the same time, Eoo[7"]- The two goals are antagonistic and, therefore, 
we are content to minimize the worst average detection delay while controlling 
the ARL to false alarm Eooir] above a prescribed level. More formally, we are 
interested in solving the minimax constrained optimization problems 

uiiJi^{T) = inf supesssupET-[(T — t)^|.7v]; subject to Eoo[r] > 7 (1-3) 

for Lorden's measure, or 

inf Jp(r) = inf supE^[T - r|T > r]; subject to Eoo[T] > 7 (1.4) 

for Pollak's measure. In both cases 7 > 1 is the prescribed minimum value of 
the ARL to false alarm. The problems (|1.3p and (II. 4p are central to sequential 
change-point detection theory, and numerous past and ongoing efforts aim to find 
the corresponding solutions for various observation models. 

Regarding existing optimality results, Lorden (1971) proved that the Cumu- 
lative Sum (CUSUM) procedure, introduced by Page (1954), asymptotically (as 
7 — > 00) solves the minimax constrained optimization problem in (jl.3p for i.i.d. 
observations before and after the change. Later, Moustakides (1986) showed that 
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CUSUM is exactly optimal for every 7 > 1 (see also Ritov (1990)). An analogous 
result for detecting a change in the drift of a Brownian motion has been indepen- 
dently established by Beibel (1996) and Shiryaev (1996). For an extension to a 
Gaussian process with independent nonhomogeneous increments see Tartakovsky 
(1995), and for a generalization to Ito processes see Moustakides (2004). 

Shiryaev (1961) introduced an alternative to the CUSUM change detection 
scheme, currently known as the Shiryaev-Roberts (SR) procedure (cf. Roberts 
(1966)), that is central here. A randomized variant of this test was proposed 
by Pollak (1985) where, instead of initializing the test from as in the original 
version, Pollak suggested a randomized initialization strategy with the initial 
point sampled from the quasi-stationary distribution of the SR statistic. We 
refer to this version as the Shiryaev-Roberts-Pollak (SRP) procedure. The gain 
obtained by this alternative initialization mechanism is significant. Pollak (1985) 
was able to demonstrate that his variant solves, asymptotically as 7 ^ 00, the 
optimization problem defined in (jl.4p within an o(l) quantity. More precisely, the 
SRP procedure has a J7p measure that differs from the (unknown) optimum by 
a quantity that tends to as 7 —> 00, even though both worst average detection 
delays tend to infinity. We refer to this asymptotic optimality as order-3, as 
opposed to order- 1, when the ratio of the two quantities tends to 1, or order-2 
when their difference is bounded. 

Despite its strong asymptotic optimality property, the SRP procedure is im- 
possible to apply in practice because there is neither an analytical nor a numerical 
method for computation of the quasi-stationary distribution required for the ini- 
tializing random variable (except in some rare cases; see, e.g., Pollak (1985) and 
Mevorach and Pollak (1991)). An important result of the present paper is a solid 
numerical technique for the computation of this distribution, making the SRP 
procedure readily available for applications. In addition, we examine alternative 
deterministic initialization strategies for the SR test. These variants, as we shall 
see in our numerical examples, are strong competitors of SRP, enjoying the same 
order-3 asymptotic optimality. As a matter of fact, in all of the examples we 
have tried, our tests (with several proposed initialization schemes) either outper- 
formed the SRP test, or performed equally well in the sense that they exhibited 
either smaller or the same J7p measure for the same ARL to false alarm. At 
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the same time, the proposed versions of the SR procedure have a fast initial re- 
sponse, providing a smaller average detection delay for changes that occur from 
the very beginning (and soon after surveillance begins) as compared to both the 
conventional SR and the SRP procedure. 

Regarding the classical SR procedure, Shiryaev (1961, 1963) considered the 
problem of detecting a change in the mean of a Brownian motion when a station- 
ary regime is in place, effected by a change possibly occurring in a distant future, 
after many false alarms have been experienced. Shiryaev proved that the SR 
procedure is exactly optimal for minimizing the expected delay in detecting such 
distant changes against a stationary background of false alarms. Recently PoUak 
and Tartakovsky (2009), motivated by this result and by the work of Feinberg 
and Shiryaev (2006), obtained a similar result for detecting a change in a general 
discrete-time model, assuming that a change occurs at a far horizon (i.e., when r 
is large) and is preceded by a stationary flow of false alarms. Specifically, the SR 
procedure was shown to be exactly optimal in minimizing, subject to the familiar 
constraint Eoo[r] > 7, the relative integral average detection delay 

instead of the worst expected conditional detection delay J-p{T) of (|1.2p . Fur- 
thermore, for a general discrete-time model, the value of J{T) has been shown 
to be equal to the limiting (as r — > 00) value of the average detection delay of 
the repeated SR detection procedure when the same stopping time is reapplied 
after each false alarm. This result was initially established by Shiryaev (1961, 
1963) for the Brownian motion model. 

Finding the appropriate version of the SR procedure that minimizes Pollak's 
i7p(T) measure is still an open problem. Answering this question is essential 
because the corresponding optimal test constitutes the missing complement of 
the CUSUM procedure for the two drastically different classes of change detection 
applications mentioned earlier. 

Concluding the literature review on the CUSUM and SR tests, Tartakovsky 
and Ivanova (1992) consider the general case of processes with independent incre- 
ments (for discrete and continuous time), providing efficient asymptotic formulas 
for the performance of the two procedures. Earlier, Pollak and Siegmund (1985) 
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carried out a similar analysis for the Brownian motion case. 

A third test is the exponentially weighted moving average (EWMA) pro- 
cedure, first proposed by Roberts (1959). Its behavior was studied in detail 
by Novikov and Ergashev (1988) and Novikov (1990) for arbitrary processes 
with independent and homogeneous increments. They show that the optimized 
EWMA procedure exhibits 23% more expected detection lag as compared to the 
CUSUM or SR procedure when detecting a change in the mean of a Gaussian 
process. These results were corroborated by Srivastava and Wu (1993) for detect- 
ing a change in the drift of a Brownian motion using an alternative technique. A 
comprehensive analysis of various EWMA schemes, and how they compete with 
CUSUM, can also be found in Lucas and Saccucci (1990). 

The main goal here is a simple numerical method for the evaluation of the 
operating characteristics of the SR test and its SRP variant. Our numerical 
technique is also used for the performance evaluation and optimization of an 
alternative version of the SR procedure in which the initializing value is deter- 
ministic instead of random (which is the case in the SRP procedure). The final 
detection procedure that comes out of this optimization (as well as its modifica- 
tions based on various initial conditions) is compared against the SRP test. In all 
numerical examples we present, the optimized deterministic initialization enjoys 
the same order-3 asymptotic optimality property as the SRP test and, more im- 
portantly, uniformly outperforms it. Of course these claims are only observations 
based on our numerical findings, but we work to support them analytically. In 
fact, a proof that the SR test with a certain deterministic initialization is exactly 
minimax (while the SRP test is not) in a particular example can be found in 
Polunchenko and Tartakovsky (2010). 

This article is organized as follows. In Section 2 we formally state the prob- 
lem and outline the SR test and its variants. In Section 3 we develop a system 
of exact integral equations on the performance metrics. In the same section wc 
propose a set of approximations to these equations that arise when we develop 
numerical solutions to the initial set of integral equations. In Section 4 we give 
numerical examples involving Gaussian and exponential models to illustrate the 
capabilities of our numerical methodology, and compare the relative performance 
of the SR test and its variants of interest. Additionally in Subsection 4.3 we show. 
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very briefly, how our computational method can be modified to suit the other 
two popular tests - CUSUM and EWMA. Finally, in Section 5 we provide a 
summary of conclusions that can be drawn from our study. 

2. The Shiryaev-Roberts Test and its Variants 

We provide a brief overview of the SR test and its randomized variant - the 
SRP test - and introduce the version with the deterministic initialization that 
we propose here as an alternative to the SRP test. 

We make the following assumptions regarding the change-point detection 
problem. Suppose a sequence {-'^n}„>i of i.i.d. random variables is observed 
sequentially. Initially the sequence is "in-control" , i.e., all the observations come 
from pdf foo{x)- At an unknown time r > 0, something happens and the sequence 
runs "out of control" by abruptly changing its statistical properties, so that from 
r + 1 on the pdf switches to /o(x) ^ foo{x)- At this point it is desired to raise an 
alarm as quickly as possible, allowing for an appropriate action to be taken. We 
recall that a sequential detection procedure is identified with a stopping time T 
that is adapted to the filtration {J^n}n>o generated by the observations. 

To define the SR procedure let = fo{Xn) / fooiXn) denote the likelihood 
ratio of the nth observation and let i?„ be the SR statistic defined as 

n n 

Rn = Y.\{^r (2-1) 

k=l j=k 

Then the original SR stopping time is the first time n that Rn attains a 
positive level i^, i.e., 

S,y = inf{n >l:Rn>iy}, with inf{0} = oo, (2.2) 

where threshold = z^^ is selected so that the false alarm constraint is satisfied 
with equality, i.e., Eoo[5jy^] = 7. It is easy to verify from (|2.ip that the SR 
statistic follows the recursion i?„ = (1 + Rn-i)in initialized with i?o = 0. 

We propose a modification by initializing the test from any value Rq = r > 0. 
Define the modified SR statistic R^ by the recursion 



i?; = (l + K_i)4, R'o = r, 



(2.3) 
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and the corresponding stopping time 



S^^ = inf{n >l:K>u} 



(2.4) 



where again u is selected so that Eoo[5^] = 7- We cah this variant SR-r. Clearly, 
threshold and initializing value r are related through the equation EqoK] = 7- 
In satisfying this equality we can either assume that is a function of r, or 
that r is a function r^, of threshold i'. For simplicity we omit subscripts. 

Our intention is to isolate a specific value for the initializing parameter r 
that will give rise to a test that competes effectively with Pollak's randomized 
SRP version. 

The SRP procedure is defined similarly to ()2.3p and (12. 4p but with Rq now 
a random variable distributed according to the quasi-stationary distribution of 
the SR statistic Rn- 



To avoid complications we assume that the likelihood ratio ii = /o(Xi)//oo(^i) 
is continuous, in which case the quasi-stationary distribution exists (cf. Harris 
(1963), Theorem III. 10.1). However, the case where ii is nonarithmetic can also 
be covered with some additional effort. 

Let the quasi-stationary density of the distribution (12. 5p be q{x). The SRP 
procedure is defined by 



where Rq ~ q{x) means that the initializing variable is random and distributed 
according to the pdf q{x). Note that q{x) = qu{x) depends on v and its support is 
[0, u). Again, the threshold v is selected so that Eoo[5^] = 7- The main drawback 
of this test has been the fact that there was no specific way to compute q{x), and 
finding q{x) has been an open problem since the first appearance of the SRP test 
in 1985. In Section 3 we give an efficient numerical answer. 

To understand the reason for the randomized initialization, we observe that 
i7p(T) is the supremum over the change time r of the sequence of conditional 
expected detection delays Et-[T — r|r > r], r > 0. According to the general 



P[i?o < x] = lim Poo[^n < x\Sl > n], x G [0,i/). 



(2.5) 




q{x) 



(2.6) 
(2.7) 
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decision theory (see, e.g., Ferguson (1967), Theorem 2.11.3) if a) we can find a 
stopping time T adapted to {J^n}n>o that is an extended Bayes and an equalizer 
rule and b) T satisfies the false alarm constraint with equality, then T solves 
the minimax constrained optimization problem defined in (jl.4p . It follows from 
Pollak (1985) that the randomized initialization according to the quasi-stationary 
distribution guarantees the equalizer property Eo[5^] = Et-[5^ — t\Su > t] for all 
r > and that threshold v = ly^ can be selected in such a way that the false 
alarm constraint is satisfied with equality. However, the SRP test was shown in 
Pollak (1985) to be only asymptotically optimal of order-3: if = z^^ is such that 
lEooi'?^] = 7, we have 

Eo[5,^l - inf supE^[r - rlT > r] = o(l) as 7 ^ oo. (2.8) 

{T:Eoo[T]>7}r>0 

The question of which test exactly optimizes i7p(T) and solves the minimax 
problem (II. 4p is still open. 

Due to ()2.8p and the fact that the SRP is an equalizer, one can conjecture 
that this test is exactly optimal. No further analysis or counterexamples were 
offered until recently to support or disprove this (see also Mei (2006) ) . We believe 
that the proposed SR-r variant can in fact provide a counterexample. As we see 
from the numerical examples in Section 4, the SR-r test, properly optimized, 
can perform uniformly better than the SRP test. Although this is only based 
on our numerical findings, it nevertheless provides a strong evidence against the 
exact optimality of the SRP procedure. A counterexample where the SRP test is 
not optimal but the proposed SR-r test is optimal can be found in Polunchenko 
and Tartakovsky (2010). While we believe that the optimal (for any given 7) 
solution is a specially designed (non-randomized) SR-r test with a varying in 
time (increasing) threshold (to guarantee constant conditional average detection 
delay), further discussion is out of the scope of this paper and will be presented 
elsewhere. 

The idea of initializing the test statistics with a value different from has 
been applied in the past by Lucas and Crosier (1982) and Lucas (1985) to 
CUSUM. The goal there was to reduce the average detection delay when ob- 
servations are affected by a change from the beginning or soon after surveillance 
begins. By means of Monte Carlo simulations, it was shown that CUSUM with 
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a positive head start exhibits the so-called fast initial response feature, permit- 
ting a more rapid response to an initial "out-of-control" situation, than does 
the conventional CUSUM (initialized from 0), at a price of a minor performance 
degradation for large values of the point of change r. However, no method for 
choosing this initial head start value has been proposed. 

2.1 Lower Bound and Asymptotic Optimality of Order-3 

To assess the quality of a detection scheme we compare the test of interest against 
a lower bound of the optimal performance. Finding such a bound turns out to 
be much easier than finding the optimal test. 

We now show that the SR-r test with initial condition Rq = r can provide 
a convenient lower bound for the optimal performance. Indeed, observe that for 
any stopping time T and any point of change r > we can write 

E.[(r-r)+] E.[(T-r)+] 
MT) >EAT-r\T>r] = = ^^^t > r] ' 

where we used the fact that since at r we are still under nominal conditions, we 
have Pt [T > t] = Pqo [T > t]. From the previous inequality we conclude that 

Jp(T)Poo[T > r] > E,[(r - r)+]. (2.9) 

Applying this inequality for r = 0, multiplying each side with r > and observing 
that Poo[^ > 0] = 1, we deduce that 

r-Jp(T) > rEo[r]. (2.10) 

Summing each side of (j2.9p over all r > and adding the corresponding sides 
of (|2.10p . we end up with the inequality 



Jp(r)ir + J^Poc[T>r] I > lrEo[T]+J2^r[{T-T)+] 

[ T=0 J [ T = Q 

or, equivalently, 

^ . rEo[T] + Er=oE.[(r-r)+] 

_ rEo[r] + Er=oJEr[(T-r)+] 
r + Eoo[r] 
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If we call the lower bound 

^ *m±i|kM^ (.„) 

r + Eoo[TJ 

and optimize each side of the previous inequality over all T that satisfy the false 
alarm constraint Eoo[r] > 7, we obtain 

inf Jp(T) > inf CviT). 

{T:E«,[T]>7} {T:E^[T]>7} 

Fortunately, the optimization of the lower bound £p(T) is possible and the opti- 
mizing stopping time is simply 5^, that is, inf|T.iE^[2-]>^} >Cp(T) = £p(5^), where 
V = \s such that EqoK] = 7- The proof of this statement for r = is given in 
Pollak and Tartakovsky (2009), and for any arbitrary positive r it can be shown 
following similar arguments. The interesting observation is that the inequality is 
true for any nonnegative value of r. 

From our previous arguments we have 

|r:Eoo[T ]>7} 

This double inequality suggests that if we are interested in verifying whether 5^ 
is asymptotically optimal of order-3, it is sufficient to show that 

hm {Jp(5:)-£p(5:)} = 0. (2.13) 
7^00 

Fix a threshold v > Q and consider the specific initializing value 

r, = arg inf {^81) - C^{S:)] (2.14) 

0<T<V 

as a candidate for initialization of the SR-r scheme. The resulting stopping time 
is now a function only of the threshold and the latter is selected so that 
SI," satisfies the false alarm constraint with equality. This uniquely defines our 
test, since both r = and u = depend only on the false alarm parameter 7. 

It turns out that the proposed initialization strategy defined in ()2.14p 
has the following property that can be used as a simpler, alternative definition. 
If we fix V and compute ^,-[5^ — t|5^ > r], r^, is the smallest r for which the 
supremum sup^>oEt-[5^ — t\SI > r] becomes equal to the steady state value 
liiRr ^00 ^ASl — t|5^ > r]. Typical forms of E,-[5^ — t|5^ > r], as functions 
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of the change time r and for different values of the initializing parameter r, 
are depicted in Fig. l2.1( a). As we can see, if r < r^, then the supremum of 
this function exceeds its steady state limit; whereas for r > r^, the supremum 
coincides with the steady state limit. We observe that 'EriSl" — ^15^" > t] 
attains the steady state value not only in the limit as r — > oo, but also for some 
finite value of r. In the same figure we can also see that the selection r = 0, 
corresponding to the classical SR test, exhibits a decreasing behavior, with the 
worst detection delay appearing at r = 0. On the other hand, the SRP test 
with threshold v has a constant performance (dashed line) which coincides with 
the steady state value. Finally, we plot the expected detection delay for r = r^, 
where is the smallest r for which Et-[5^ — t|5^ > r] becomes an increasing 
function of r. 




Change time r Initializing variable r 

(a) (b) 

Figure 2.1: (a) Typical form of expected detection delay as a function of change-time 
r for various initialization strategies and (b) ARL to false alarm as a function of the 
initializing parameter r. 



Therefore, we conclude that, for any given threshold u, the proposed initial- 
izing parameter ri, can be alternatively defined as 

= inf |r > : - A^l >t]< lim - t\SI > r]| . (2.15) 

L r— >oo J 

The advantage here is that this does not involve the computation of the lower 
bound C-p{Sl). 

In Fig. l2.1l fb) we plot the ARL to false alarm EqoK] as a function of the 
initializing parameter r. For the same threshold, the proposed SR-r^, test exhibits 
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a larger ARL to false alarm than the one obtained by the SRP test. Since both 
tests have the same worst conditional average detection delay (i.e., the same 
Pollak measure), this suggests that initialization with r = r^ is preferable to the 
SRP. Note that all values of r inside the half-toned strip in Fig. l2.1f b) correspond 
to tests that perform better than the SRP procedure. Furthermore, since the 
corresponding values of r are larger than r,^, their worst detection delay is equal 
to the steady state value and thus equal to the SRP performance. 

As a last remark, recall that corresponds to the smallest r for which the 
conditional average detection delay becomes an increasing function of r. Thus 
this initialization strategy exhibits a strong fast initial response, responding much 
faster to changes that take place in earlier than later stages, complete opposite 
to the classical SR test that starts from r = and prefers large change-times r. 

3. Proposed Methodology 

We derive here exact integral equations, as well as relevant recursive formulas, 
for the performance metrics of the SR-r and SRP tests. The exact formulation 
is then followed by a set of approximations leading to classical linear problems 
from Linear Algebra that can be solved numerically, and thus provide answers 
to long standing performance evaluation problems. 

3.1 Integral Equations for the Operating Characteristics of the SR-r 



Fix r, u with r £ [0, v) and define (j)i{r) = Ej[5^], where i = 0, oo. The function 
0oo(r) is the ARL to false alarm and (po{r) is the average detection delay when 
the change takes place before surveillance begins. Since the statistic i?^ obeys the 
recursion = {1 + R'^_i)£n (cf. (j2.3p ). it is clear that {Rn}n>o is a homogeneous 
Markov process and, therefore, one has 



where -Fj(-) is the cdf of the likelihood ratio ii under the Pj measure, by substi- 
tuting this equality into the previous one we obtain the final integral form for 



Test 



= r 



where, hereafter, 1_4 stands for the indicator of a set A. Since 
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the functions of interest 

Mr) = '^ + (t>i{^)-^F'i(^j^^ dx, i = 0,oo. (3.1) 

The next important performance metric is the conditional expected detection 
delay 

E,[5:-r|5:>r] = %M|^^, r = 0,l,2.... (3.2) 

We consider the numerator and the denominator separately and propose recursive 
formulas for their computation. For r > 1, let Srif) = lEr[(5^— t)"'"] and Prir) = 
PooK > ''"]• -Du^ to Markov nature of {i?5^}„>o, for r > 1, we have recursions 
for the sequences of functions: 

= L ^^-'^''^ ^""-^ {tt^) 

Pr{r) = Eoo[pr-l{Rl)l{Rl<u}\Ro = r] 

where 6o{r) = (poi^r) and po{r) = 1. The sequences of functions {ST{r)},{pr{r)} 
can be computed using these recursive formulas; they involve a repetitive appli- 
cation of the same linear transformation on an initial function and 

EAS: - t\S: >r] = ^, JpiS:) = sup ^^^''^ 



Prir) " r>0 Prir) 

Remark. It is of interest to evaluate the local false alarm probabilities Poo [•S^ < 
k + m\Sl, > k], k = 0, 1, ... , inside a fixed "window" of size m (m > 1; for 
m = 1 we obtain the instantaneous false alarm probability). In particular, the 
supremum local false alarm probability sup^>oPooK ^ k + m\Sl, > k] can serve 
as an alternative measure of the false alarm rate in place of the ARL to false 
alarm (see Tartakovsky (2005, 2009) for a more detailed discussion). Since 

r^[s:<k+m\s:>k] = l^\^^^^±^ = i ^-K>^+-] 



we obtain that PooK < ^ + nT-lSl, > A;] = 1 — Pk+m{f)/Pk{f)^ where po(^) = 1 
and Pj{r), j = 1,2 . . . are given in (|3.4p . 
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3.2 Integral Equations for the Operating Characteristics of the SRP 
Test 

We now derive equations for the operating characteristics of the randomized SRP 
scheme. First, the quasi-stationary density q{x) satisfies the integral equation 

Amax9(2;) = / q{r)-^Foo ( dr (3.5) 
Jq dx \l + rj 

(cf. Pollak (1985)), where A 

max is the leading eigenvalue of the linear integral 
operator induced by the kernel 

d I X \ 
Kocix, r) = — Foo ( ) , X, r € [0, z^), 

and, consequently, q{x) is the corresponding (left) eigenfunction. Since q{x) is a 
probability density with support [0,u), it also satisfies the constraint 

q{x) dx = 1. (3.6) 



Equations (j3.5p and (j3.6p are sufficient to uniquely define Amax and q{x), while 
their existence is guaranteed by Harris (1963), Theorem III. 10.1. Furthermore, 
integrating both sides of ()3.5p with respect to x over the interval [0, u), using ()3.6p 
and the fact that -Foo(O) = 0, we conclude that 

< Amax = ^ q{r)F^ (iT^) ^ / '^^'^^ "^"^ " ^' 

Note that Foo{i^/{l + r)) < 1, but there is an interval for r where this inequality 
is strict thus implying the strict inequality in the previous relation. Indeed, by 
assuming that the pdfs before and after the change are different, we have that 
-Poo(l) < 1- Using the continuity of -^00(2;) (as a result of the assumption that 
li is continuous) we have that Foo{x) < 1 for x < 1 and sufficiently close to 1. 
Now note that, for sufficiently large i' and r sufficiently close to z^, we assess that 
1^/(1 + r) < 1 and that it is sufficiently close to 1, therefore Foo(i//(l + r)) < 1. 
Thus the leading eigenvalue Amax is nonnegative and strictly bounded by 1, and 
the same is true for Fo{x/{l + r)). 

Assuming that q{x) is available through a solution of (j3.5p . we proceed with 
the computation of the performance of 5^. Since S$ is an equalizer rule, c7p(5^) = 
]Eo[5^], while the ARL to false alarm becomes Eooi^t^]. In the case of the SRP 
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test, averaging is with respect not only to the observation statistics, but also to 
the distribution of the initializing point. The expected values ]Eoo[5^] and ]Eo[52] 
are easy to compute when the quasi-stationary pdf q{x) and the two functions 
(j)iir), i = 0,00, are available. Indeed, 

Ei[S^,]= [ Ei[S^\Rl = x]q(x)dx = [ Ei[S^]q{x)dx 

Jo Jo (3 7) 

(pi(x)q(x)dx. 





3.3 Lower Bound Computation 

Using the definition of Srir) in (j3.3p and the fact that the leading eigenvalue 
of the linear transformation that updates the sequence of functions {Sr{r)} is 
'^max) we conclude that d-rir) = 0{X'^^^). Since < Amax < !> it follows that the 
series in the numerator of the lower bound YlT=o^'r[i'^u ~ '^)~^] ~ Yl'^=o is 
absolutely summable. Consequently, it is easy to verify that ^p{r) = X^!^o^'^(^) 
is the solution of the integral equation 

Using (j2.1ip and the above notation, we obtain 



where (piir), i = 0,oo are given by (j3.ip and ^(r) by (j3.8p . 
3.4 Numerical Solutions 

Observe first that (|3.ip and (|3.5p can be written in the form 

u{r) — a / K{x,r)u{x) dx = v{r), (3.10) 



Jo 

where u{x) is an unknown function, a ^ and x,r £ [0, u]. In particular, (j3.ip 
can be obtained from (j3.10p by letting u{x) = (j)i{x), K{x,r) = -^Fi{x/{1 + r)), 
i = 0,00, v{r) = 1, and a = 1. Replacing integration over x in (jS.lOp with 
integration over r, and u{x) with u{r) = q{r) under the integral and choosing 
v{r) = 0, K{x,r) = ^Foo(x/(l + r)), and a = 1/Amax, yields (|3.5I) . We assume 
sufficient smoothness of our functions so that the interval [0, z^) can be extended 
to [0,zv]. 
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An equation analogous to (j3.10p occurs in a wide variety of physical applica- 
tions and is known as the Predholm equation of the second kind; see, e.g., Petro- 
vskii (1957) and Kress (1989). The fundamental result concerning the existence 
and uniqueness of solutions of such equations is that, given certain regularity 
conditions on the kernel, these equations have unique solutions provided 1/a is 
not a proper number or an eigenvalue of the linear integral operator associated 
with the kernel K[x, y). As we have seen, a = 1 is not an eigenvalue of either of 
the operators induced by -^Fi{x/{1 + r)), i = 0, oo. 

Various numerical schemes for solving (jS.lOp are developed in Kantorovich 
and Krylov (1958), Petrovskii (1957), and Atkinson and Han (2001). Commonly, 
one replaces the function /(r) = K{x,r)u{x)dx in (|3.10p by a vector / = 
[/(ro), /(n), . . . , /(fat)]*, where = ro < ri < ••• < tat = u constitutes a 
sampling of the interval [0, z^]. A similar sampling is applied to the function 
u{x) producing the vector u = [u{xq),u{xi), . . . ,u{x]\f)Y . The integral is then 
evaluated using some numerical integration technique, leading to a (right) matrix- 
vector multiplication that replaces the integral, 

f{r)=[ K{x,r)u{x)dx ^ f = Ku, (3.11) 
Jo 

where is a matrix that depends on the numerical integration method and the 
sampling points {rj, {xj and / = [/(ro), /(ri), . . . , /(rAr)]*, with /(r) denoting 
the approximation to /(r) as a result of evaluating the integral numerically. The 
same matrix K used for the numerical evaluation of the integral in (13. lip can 
also be used to evaluate the conjugate integral via 

f{x)= [ K{x,r)u{r)dr ^ f = u^K. (3.12) 
Jo 

To find the matrix K we need to use numerical integration. We employ the 
simplest numerical technique to demonstrate the main idea; if one adopts more 
powerful numerical integration methods, the results will be of higher accuracy. 

Consider an integral of the form J^^ z{x)F' (x) dx, where F'{x) is the deriva- 
tive of F{x) and let a = < xi < • • • < xat = 6 be a sampling of the interval 



18 G.V. MOUSTAKIDES, A.S. POLUNCHENKO AND A.G. TARTAKOVSKY 

[a,b]. Then 



.6 1 ^ 

/ z(x)F'(x)<ix« -^[F(x,-)-F(x,_i)][z(x,) + 

^ .7 = 1 



-[F{xi) - F{xo)]z{xo) 



N-l 



(3.13) 



+ ^[F{xn) - F{xn-i)]z{xn). 

For (j3.ip and (j3.5p ). as well as the computations in (j3.3p and (j3.4p that 
involve the kernels Ki{x,r) = ^-^iCxf^)) ^ = 0, oo, sample the interval [0, i^] 
canonically at Xj = Vj = jc, j = 0, . . . , N with c = ly/N, and apply (|3.1ip . (|3.12p 
and (j3.13p to obtain 



/' 

Jo 



" d_ 

dx 



F 



X 



1 + r 

X 



dx \1 + r 



u{x) dx =^ MiU, 
u{r) dr =^> u^Ni. 



(3.14) 
(3.15) 



The matrices iWj, Ni are of size (A^ + 1) x (A^ + 1) with elements 



{Mi 



<k,m, 



0.5Fi 
n 5/?. littlk 



0.5K 



l+mc 

- 0.5F, 



l+mc 



(fc-l)c ^ 
* I l+mc 
{N-l)c^ 



for /c = 0, 

for iV > A; > 1, 



(3.16) 



0.5F,- for A: = N, 

' ' l+mc / ' 



where in the first line of (j3.16p we used the fact that -Fj(O) = 0; and 



(AT, 



i)m,k 



0.5Fi 

0.5K; 



O.SFi {mc)-0.5Fi (^f 



mc 



l+{fc+l)c 
mc 



0.5K 



0.5Fi 



\^l+(fc-l)c 



for k = 0, 
iov N > k > I, 
for k = N. 



(3.17) 



l+(Af-l)c^ 

Using ([3Ti]) . (f3T]) and (fS^Sj) are reduced to 

= J + Mi4>i, i = 0, oo, 

= 00 + MooV', 

where 0^ = [(^^(0), <^j(c), • • • ,0i(j^)]*, = [^^(0), V5(c), • • • ,^(i^)]*, with (^^(x), 
'(/'(x) denoting the approximation to (j)i{x) and ip{x), respectively, and J = 



(3.18) 
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[1 • • • 1]*. Solving the linear system of equations in ()3.18p yields the required 
approximation for the functions (j)i{r) and ^(r). 

The recursions in (13. Sh and (13. 4p . using again (|3.14p . can be approximated 

by 

K = M^K^i, ~So = 4, (3.19) 
= MooP,_i, Po = J, (3.20) 

where Sr = [6r{0),6r{c),5r{2c),- ■ ■ ,5r{i^)Y, c = v/N, and 5t-{x) denotes the 
approximation to 6r{x). A similar definition applies to p^. 

We now turn to (|3.5p . Since this involves conjugate integration, we need to 
apply the approximation given in (j3.15p . which leads to 

AmaxQ* = q'N^. (3.21) 

This suggests that (AmaxjQ) is a left eigenvalue-eigenvector pair for the matrix 
Noo with Amax being the leading eigenvalue of N^o- Of course q is not unique 
unless we use the constraint (j3.6p . Applying (|3.13p , this constraint is transformed 
to 

c[0.5, I,-- - ,1, 0.5]q = 1, (3.22) 

where c = v/N. Since the matrix has positive elements, see ()3.17p . its 
leading eigenvalue Amax and the corresponding left and right eigenvectors (conse- 
quently also q) are necessarily nonnegative (see, e.g., Horn and Johnson (1990)). 
Following similar arguments as in Subsection 3.2 we can show that < Amax < 1- 
For computing the performance of the SRP procedure from (j3.7p we have 

%[Sl] « c[0.5, 1, • • • , 1, O.5](0, o q), (3.23) 

where, if x, y are vectors of the same length, xoy denotes the vector that results 
from the element-by-element multiplication. 

Now p.lSp - ()3.23p can be used to obtain numerical solutions for our perfor- 
mance evaluation problem. When v is large and/or sampling is fine, solving (j3.18p 
and ()3.2ip can be particularly memory and time demanding. For such large-sized 
problems it is preferable to apply iterative solution techniques that avoid storage 
of the matrices Mi,Ni and, for Eq. (j3.18p . to use simple pre-conditioning tech- 
niques to speed up convergence (see, e.g., Quarteroni, Sacco, and Saleri (2000)). 
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Approximation accuracy is, of course, directly related to A^. If is suffi- 
ciently large, the numerical solution is close to exact, see Kantorovich and Krylov 
(1958) or Atkinson and Han (2001). More details on computation of efficient up- 
per bounds for the numerical errors are reported in Polunchenko (2009). 

In the next section, we apply these ideas in specific examples in order to 
compare the relative performance of different initialization strategies. 

4. Numerical Examples 

Apart from the initialization strategies introduced in Subsection 2.1, we also ex- 
amine the case of SR-/i with r = fi = xq[x)dx, where we initialize the test with 
the mean of the quasi-stationary distribution. Recall that SR-r^ corresponds to 
the smallest value of r for which Et-[5^ — t|5^ > r] becomes increasing with 
respect to r (for fast initial response) and SR-r,y is the initialization obtained 
by solving ()2.14p or p.lSp . We did a comparative study of the change detection 
procedures for two models - Gaussian and Exponential. 

4.1 Gaussian Example 

Consider a Gaussian example of detecting a change in the mean value where 
observations are i.i.d. M{0, 1) pre-change and i.i.d. M{9, 1) post-change: 



We performed extensive numerical computations for various parameter values, 
but present only sample results for 9 = 0.1, corresponding to a relatively small 
change that is not easily detectable. We considered Eoo[5^] = 10^, 10^, as mod- 
erate and low false alarm rates, corresponding to detection threshold v in the 
range 1000 it 10% and 10000 it 10%, respectively. The integration interval [0, i^] 
was sampled at A^ = 10^ and A^ = 10^ equidistant points. We believe that such 
sampling is sufficiently fine since the results of Monte Carlo experiments for the 
conventional SR procedure (with 10^ replications) matched our numerical results 
within 0.5%. 

It is important to have a fairly accurate initial guess in order to obtain 
a pilot estimate of Eoo[5^] in searching for appropriate threshold values in a 
relatively narrow interval. To this end, the approximation Eoo[5^] = i^/w — r is 
used, where the constant w G (0,1) (related to the "overshoot") is the subject 




and 





Figure 4.1: Average detection delay Er[5^ — t|5^ > t] for the different initialization 
strategies as a function of t for 9 = 0.1 and for ARL to false alarm: (a) 7 = 10'^, 
(b) 7 = 10^. 



of renewal theory and can be computed numerically. This approximation can be 
obtained by noticing that — n — r is a PcxD-martingale with zero expectation. 
Consequently, by the optional sampling theorem, we have KoolR^r — S^ — r] = 0. 
Hence EqoK] = Eooi-R^r] — r and, since K^r is the first excess over u, renewal 
theory can be applied to the "overshoot" log(-R^r) — logi^. This approximation 
was first derived for r = in Pollak (1987), and its generalization for any r £ 
(0, u) is straightforward. For the SRP procedure the value of r should be replaced 
by /i = Eoo[-Ro], the mean of the quasi-stationary distribution. 

Fig. 14.1( a) shows the family of curves E,-[5^ — t|5^ > r] versus r for all 
initialization procedures in question when 9 = 0.1 and ARL to false alarm 7 = 
10"^. Fig. l4.1l fb) depicts the same curves for 7 = 10^. In order to have a more 
precise idea of the relative performance difference of the competing schemes, in 
TableHTT] we list the numerical values obtained by our computational method 
for characteristic values r and the ARL to false alarm of 10^. Table lT2] depicts 
the thresholds and the corresponding values of the initializing parameter r that 
assure the desired values of the ARL to false alarm for each initialization strategy. 

As we can see, the SRP procedure maintains constant average detection 
delay as expected. The SR-r* test has the fastest initial response (for immediate 
and early changes), but the worst minimax behavior. The SR-r,y procedure is 
uniformly better than the SRP test. Even though the difference is not dramatic 
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Table 4.1: Average detection delay — t|5^ > t] versus change point r for the 

ARL to false alarm 7 = 10^ and 9 = O.l. 



Tcst\r 





50 


100 


200 


400 


600 


800 


1000 


SR 


298.5 


258.3 


230.2 


197.7 


182.9 


181.5 


181.4 


181.4 


SR-r,y 


202.8 


195.9 


196.4 


200.1 


202.5 


202.8 


202.8 


202.8 


SR-n 


174.9 


179.9 


191.6 


205.6 


213.1 


214.1 


214.2 


214.3 




194.0 


190.7 


194.6 


201.6 


205.6 


206.0 


206.1 


206.1 


SRP 


206.1 



Table 4.2: Detection thresholds and initializing parameters resulting in the ARL to false 
alarm 7 = lO^, 10^ for 9 = 0.1. 



7 


10^ 


10^ 


Test 


V 


r 


V 


r 


SR 


944.0 





9435.0 





SR-r^ 


1142.0 


210.8 


9775.0 


361.2 


SR-r, 


1258.0 


333.2 


9792.0 


380.4 


SR-^ 


1174.0 


244.4 


9945.0 


540.9 


SRP 


1174.0 


random 


9945.0 


random 



it is nonetheless visible here. It is interesting to note that the SR-// detection 
procedure has an intermediate performance between SR-rj, and SR-r^, namely 
sufficiently fast initial response and the same minimax performance as the SRP 
test attained at steady state. 

Regarding the conventional SR test (with r = 0) note that it outperforms 
all competing schemes including SRP for sufficiently large change-time r. This 
is expected since, as can be seen from Fig. l2.1( b). when all tests have the same 
threshold the SR test has the largest ARL to false alarm and the same steady- 
state value for the expected detection delay. In the other tests, in order to attain 
the same ARL to false alarm value as SR, thresholds should be increased. This 
will result in an increase in the expected detection delay and, in particular, the 
corresponding steady state value. Consequently, the expected delay of SR, due to 
its monotone behavior, attains smaller values than the other tests for sufficiently 
large change-time r. 

Fig. l4.2l (a) and (b) depict the supremum average detection delay ^p(5^) as 
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ARL to false alarm 7 ARL to false alarm 7 

(a) (b) 



Figure 4.2: Worst average detection delay of competing initialization strategies and 
corresponding lower bound as a function of the ARL to false alarm 7 for = 0.1, 
Gaussian case. 

a function of the ARL to false alarm Eoo[5^] for the initialization strategies of 
interest, along with the lower bound /^p(5^''). We can see that the SR-rj^ test 
uniformly outperforms all its rivals. We also observe that initializing the SR test 
deterministically with the mean of the quasi-stationary distribution results in a 
performance that is indistinguishable from that of SRP. Finally, we can see that 
the SR-r^ procedure is inferior to SRP, but we recall that this version of the SR 
test has the best performance in terms of fast initial response. Summarizing, the 
best minimax performance is delivered by the SR-r,y test. This performance is 
also very close to the lower bound £p(5^''), suggesting that the unknown optimal 
test can offer only insignificant improvement over SR-r,y. 

4.2 Exponential Example 

Consider now the case where observations are independent, originally having an 
Exponential(l) distribution, changing at an unknown time r to Exponential(0), 

/oo(2;) =e-"l|,>0}, /o(a:) = e-'e-fl|,>o}, ^ > 0. (4.1) 

Fig. l4.3l (a) and (b) depict the supremum average detection delay J-p{Sl) 
versus the ARL to false alarm for the SRP and SR-r,y detection procedures, 
along with the lower bound C-p{Sl/') when 9 = 1.1. As in the Gaussian example, 
the operating characteristic of the SR-r,y procedure is uniformly better than that 
of the SRP procedure, with the worst average detection delay J'p{Sl'') being very 
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150| < < 1 < 1 210 




ARL to false alarm 7 ARL to false alarm 7 

(a) (b) 

Figure 4.3: Worst average detection delay of SRP and SR-ri, and corresponding lower 
bound as a function of the ARL to false alarm 7 for 9 = 1.1, exponential case. 

close to the lower bound of the minimax risk. Again we observe that there is 
very little margin for improvement over the proposed detection procedure SR-rj,. 

4.3 Extensions for CUSUM and EWMA 

Extending our results to cover the case of CUSUM and EWMA presents no 
special difficulty. Consider first the CUSUM procedure defined by the CUSUM 
statistic and the corresponding stopping time as 

V; = max{l, V:_,}en, V^^ = r; = inf{n > 1 : K' > 

where < r < u. Note that the classical CUSUM is initialized with r = 1 (cf., 
e.g., Moustakides (1986)). Here we consider the more general case suggested 
by Lucas and Crosier (1982), which leads to fast initial response. When we 
compare (I2.3p and (12. 4p with the previous two formulas, we find a difference 
only in the term (1 + R'^_i), which is now replaced by max{l, V^_]^}. The same 
difference is encountered in the integral equations that specify the two functions 
(j)i{r) = Ki[TJ'], i = 0,00, as well as the recursions that compute the numerator 
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(j)i{r) = 1+1 4'i{x)-^Fi [ ^- — - ] dx, i = 0,oo, 



X 


max{l, 


r} 


X 




max{l, 


r} 


X 




max{l. 


r} 



and denominator of IEt['?7 — > ''"]• Specifically. 

d_ 

dx 

^rir) = [ 6r-i{x) -^Foo ( ^- — 7 ) dx, 6o{r) = (poir), 

Jq ox ymaxjl, rjy 

Pr{r) = [ pr-i{x) -^Foo ( Y^ — T ) P^^^^ = ^• 

Jo ox \max{l, r| / 

A randomized version with Vq following the quasi-stationary distribution P[Vo < 
x] = lim„^oo ^[Vn ^ ^\%} > ^] is also possible to define, with pdf satisfying 

In fact randomization is not necessary since the conventional CUSUM test with 
r = 1 is exactly optimal in the sense of Lorden, (jl.ip . and the randomized 
CUSUM is always inferior to the SRP test in the sense of Pollak, ()1.2p . 

Approximations that produce numerical solutions to these equations can be 
found in the same way. Note that Dragalin (1994) used a slightly different, 
but very precise, numerical technique for the computation of the ARL to false 
alarm Eooi?^^] and the average detection delay Eo[T^] of the standard CUSUM 
for the Gaussian distribution. Comparison of results obtained by our numerical 
technique with those obtained by Dragalin shows that the two approximations 
are close, suggesting that our simple numerical method is of sufficiently high 
precision. 

Finally, in order to make a similar extension for EWMA, here 

where is the EWMA statistic, A/"^^ corresponding (double sided) 

stopping time, < a < 1 is a forgetting factor, and < i^i < 1 < 1^2 ai'e the 
thresholds (the case 1^1 = corresponds to the one-sided EWMA procedure). 
Note that the EWMA statistic is usually written in a form involving the log- 
likelihood ratio log(£„). This conventional form can be recovered by simply taking 
the logarithms, but we prefer the exponential version above, since it allows for the 
derivation of integral equations for a variety of performance metrics. Observing 
that the difference with the SR case is that the term (1 + Rn-i) replaced by 
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{D^_^)^, yields the equations that define the operating characteristics 

f"^ d / X \ 
^i{r) = 1 + / '^»(^)^-^» (^j ^ = 0>oo, 

6^{r) = / 6r-iix) -Q^F^ [^j dx, 6Q{r) = (po{r), 

Prir) = J pT-i{x) ^Foo [— ) dx, po{r) = 1, 
f"^ d / X \ 

Amaxg(x) = J^^ 9(r)^Foo ( — J dr, 

where the last expression corresponds to the pdf of the quasi-stationary distribu- 
tion defined as P[£>o < x] = lim„^oo IPoo [-DJ^ < x\Nl^^^^ > n] when a randomized 
EWMA is being constructed. To our knowledge, no randomized EWMA scheme 
has been previously considered. 

Producing numerical approximations is again straightforward. Note that 
Robinson and Ho (1978) proposed a different approach to obtaining numerical 
approximations for the performance of a somewhat different EWMA procedure. 
We believe that our approach is advantageous because it allows not only for 
the evaluation of the ARL to false alarm EooiA/"^^ j^^] ™d the ARL to detection 
Eo[A/'^^^i^2], but also for the optimization of the initializing parameter r. Fur- 
thermore, our method can be used to find the change time r that produces the 
worst conditional expected detection delay. As opposed to the standard SR and 
CUSUM, in this test it is expected that this worst performance appears at a time 
r > 0. 



5. Conclusion 

For the problem of quickest change detection with known pre- and post-change 
distributions, we proposed a simple modification of the SR procedure, called the 
SR-r test, that starts from a deterministic (fixed) point r > 0. This procedure 
represents a family of sequential tests, as the initializing point r can take any 
value in the interval [0, i/), with v denoting the threshold. Our main contribu- 
tion is the development of integral equations for the major performance metrics 
of the SR-r test and the corresponding numerical techniques for solving these 
equations. Additionally, we give a method for the numerical computation of the 
quasi-stationary distribution of the SR statistic. This allows the practical imple- 
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mentation of the randomized SRP procedure, introduced by Pollak (1985), that 
is known to enjoy strong asymptotic optimahty properties. Using our numerical 
methodology we can compute the operating characteristics of the conventional 
SR procedure (that starts from 0) and of several interesting SR-r variants cor- 
responding to specific choices of the initializing parameter r, and can compare 
their performance against the SRP test. 

The numerical results, obtained for the Gaussian and Exponential exam- 
ples, indicate that a specially designed SR-r test can uniformly (for all points of 
change) outperform the SRP procedure. Even though the difference is not dra- 
matic, this observation constitutes strong evidence against the exact optimality 
property of the SRP procedure. As Polunchcnko and Tartakovsky (2010) prove, 
using the integral equations presented in Subsections 3.1 and 3.2, the SRP proce- 
dure is indeed not strictly optimal, while the SR-r procedure may be optimal in 
certain examples. Our numerical analysis also shows that by slightly sacrificing 
performance in the worst case detection delay sense, it is possible to design SR-r 
tests that exhibit fast initial response (i.e., guarantee faster detection of changes 
that occur at early stages). 
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